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Detailed Fermi-surface structures are essential to describe the upper critical field H C 2 in type- 
II superconductors, as first noticed by Hohenberg and Werthamer [Phys. Rev. 153, 493 (1967)] 
and shown explicitly by Butler for high-purity cubic Niobium [Phys. Rev. Lett. 44, 1516 (1980)]. 
We derive an H C 2 equation for classic type-II superconductors which is applicable to systems with 
anisotropic Fermi surfaces and/or energy gaps under arbitrary field directions. It can be solved 
efficiently by using Fermi surfaces from ab initio electronic-structure calculations. Thus, it is ex- 
pected to enhance our quantitative understanding on H C 2- Based on the formalism, we calculate 
Hc2 curves for Fermi surfaces of a three-dimensional tight-binding model with cubic symmetry, an 
isotropic gap, and no impurity scatterings. It is found that, as the Fermi surface approaches to the 
Brillouin zone boundary, the reduced critical field h*(T/T c ), which is normalized by the initial slope 
at T c , is enhanced significantly over the curve for the spherical Fermi surface with a marked upward 
curvature. Thus, the Fermi-surface anisotropy can be a main source of the upward curvature in H C 2 
near T c . 

PACS numbers: 74.25.Op, 71.18.+y 



I. INTRODUCTION 

The upper critical field H C 2 is one of the most fun- 
damental quantities in type-II superconductors. Af- 
ter the pioneering work by Abrikosov 1 based on the 
Ginzburg-Landau (GL) equations, 2 theoretical efforts 
have been made for its quantitative description at all 



temperatures 



3-41 



However, we still have a limited suc- 



cess when compared with those for the electronic struc- 
tures in the normal state. 42 The purpose of the present 
paper is to provide a theoretical framework which enables 
us ab initio calculations of as accurate as electronic- 
structure calculations in the normal state. 

Necessary ingredients to be included are (i) nonlocal ef- 
fects effective at low temperatures; (ii) impurity scatter- 
ing; (hi) Fermi-surface anisotropy; (iv) strong-coupling 
effects; (v) gap anisotropy; (vi) mixing of higher Lan- 
dau levels in the spatial dependence of the pair poten- 
tial; (vii) Landau-level quantization in the quasiparti- 
cle energy; 13 ' 32 ' 35-38 (viii) fluctuations beyond the mean- 
field theory. 43 We here derive an H C 2 equation which 
is numerically tractable, including all the effects except 
(vii) and (viii). 

An H C 2 equation considering the effects (i) and (ii) 
was obtained by Helfand and Werthamer. 6 It was ex- 
tended by Hohenberg and Werthamer 9 to take the Fermi- 
surface anisotropy (iii) into account. Equations with the 
strong-coupling effects (iv) were derived by Eilenberger 
and Ambegaokar 11 using Matsubara frequencies and by 
Werthamer and McMillan 12 on the real energy axis, 
which arc equivalent to one another. Schossmann and 
Schachinger 27 later incorporated Pauli paramagnetism 
into the strong-coupling equation. Although an equa- 



tion including (i)-(iv) was presented by Langmann, 33 it 
is still rather complicated for carrying out an actual nu- 
merical computation. On the other hand, Rieck and 
Scharnberg 30 presented an efficient H C 2 equation where 
the effects (i)-(iii) and (vi) were taken into account, and 
also (v) in the special case of the clean limit. See also 
the work by Rieck, Scharnberg, and Schopohl 31 where 
the strong-coupling effects (v) have also been consid- 
ered. Our study can be regarded as a direct extension 
of the Rieck-Scharnberg equation 30 to incorporate (i)- 
(iv) simultaneously. To this end, we adopt a slightly 
different and (probably) more convenient procedure of 
using creation and annihilation operators. We will pro- 
ceed with clarifying the connections with the Rieck- 
Scharnberg equation as explicitly as possible. 

The remarkable success of the simplified Bardeen- 
Cooper-Schrieffer (BCS) theory 44,45 tells us that detailed 
electronic structures are rather irrelevant to the proper- 
ties of classic superconductors at H = 0. However, this 
is not the case for the properties of type-II supercon- 
ductors in finite magnetic fields, especially in the clean 
limit, as first recognized by Hohenberg and Werthamer. 9 
Their effort to include the Fermi-surface anisotropy in 
the H C 2 equation was motivated by the fact that the 
Helfand- Werthamer theory 6 using the spherical Fermi 
surface shows neither qualitative nor quantitative agree- 
ments with experiments on clean type-II superconductors 
like Nb 46 ~ 48 and V. 49 Indeed, angular variation in H C 2 by 
10% was observed at low temperatures in high-quality 
Nb 46 ' 50 ' 51 and V 50,51 with cubic symmetry. 52 Also, the 
reduced critical field 



h*(t) 



H c2 {t) 



-dH c2 (t)/dt\ t=1 



(t = T/T c ) , 



(1) 
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calculated by Helfand and Werthamer 6 has h*(0) = 0.727 
in the clean limit, whereas a later experiment on high- 
purity Nb shows (h*(0)) = 1.06 for the average over 
field directions. 51 Hohcnbcrg and Werthamer 9 carried 
out a perturbation expansion for cubic materials with re- 
spect to the nonlocal correction where the Fermi-surface 
anisotropy enters. They could thereby provide a quali- 
tative understanding of the H c2 anisotropy and the en- 
hancement of (h* (t)} observed in Nb. They also derived 
an expression for (h* (0)) applicable to anisotropic Fermi 
surfaces. It was later used by Mattheiss 14 to estimate 
(h*(0)) = 0.989 for Nb based on his detailed electronic- 
structure calculation. The strong dependence of h* (t) in 
the clean limit on detailed Fermi-surface structures can 
also be seen clearly in the numerical results from a model 
calculation by Ricck and Scharnberg, 30 and from the dif- 
ference h*(0) = 0.727 and 0.591 between spherical and 
cylindrical Fermi surfaces, respectively. 41 

On the other hand, it was shown by Werthamer and 
McMillan 12 that the strong-coupling effects change h* (t) 
by only <2% for the spherical Fermi surface and cannot 
be the main reason for the enhancement of h* (0) in Nb. 

The most complete calculation including the effects 
(i)-(iv) was performed on pure Nb by Butler. 22 ' 23 He 
solved the strong-coupling equation by Eilenbcrger and 
Ambegaokar, 11 taking full account of the Fermi-surface 
structure and the phonon spectra from his electronic- 
structure calculations. He could thereby obtain an excel- 
lent agreement with experiments by Williamson 50 with 
(h* (0)) = 0.96 and by Kerchner et a/. 53 However, a later 
experiment by Sauerzopf et al. on a high- purity Nb 
shows a larger value (h*(0)) = 1.06, thereby suggesting 
that there may be some factors missing in Butler's cal- 
culation. 

Theoretical considerations on the effects (v) and (vi) 
started much later. It was Takanaka 18 and Teichler 19 ' 20 
who first included gap anisotropy (v) in the H C 2 equation. 
They both considered the nonlocal effect perturbatively 
adopting a separable pair potential. Takanaka studied 
H C 2 anisotropy observed in uniaxial crystals, whereas Te- 
ichler applied his theory to the H C 2 anisotropy in cubic 
Nb. This approach by Teichler was extended by Proham- 
mer and Schachingcr 28 to anisotropic polycrystals and 
used by Weber et al. 54 to analyze anisotropy effects in 
Nb. 

The mixing of higher Landau levels (vi) was first 
considered by Takanaka and Nagashima 15 in extending 
the Hohenberg- Werthamer theory for cubic materials 9 
to higher orders in the nonlocal correction. It was 
also taken into account by Takanaka 18 in the above- 
mentioned work, by Youngner and Klemm 24 in their per- 
turbation expansion with respect to the nonlocal cor- 
rections, by Scharnberg and Klemm 25 in studying H C 2 
for p-wave superconductors, by Rieck and Scharnberg 30 
for superconductors with nearly cylindrical model Fermi 
surfaces, and by Prohammer and Carbotte 29 for ci-wave 
superconductors. See also a recent work by Miranovic, 
Machida, and Kogan on MgB 2 . 39 Although it plays an 



important role in the presence of gap anisotropy, 25 ' 29 this 
mixing was not considered by Teichler. 19 ' 20 

Now, one may be convinced that calculations includ- 
ing (i)-(vi) arc still absent. Especially, many of the 
theoretical efforts have been focused only on the spe- 
cial case of cubic materials. 9 ' 15 ' 19 ' 20,22 ' 23 For example, 
a detailed theory is still absent for the large positive 
(upward) curvature observed in H c2 (T < T c ) of lay- 
ered superconductors, 55 ' 56 except a qualitative descrip- 
tion by Takanaka 18 and Dalrymple and Prober. 57 Based 
on these observations, we here derive an H c2 equation 
which is numerically tractable for arbitrary crystal struc- 
tures and field directions by using Fermi surfaces from 
ab initio electronic-structure calculations. This kind of 
calculations has been performed only for Nb by Butler 
so far. 22 ' 23 Making such calculations possible for other 
materials is expected to enhance our quantitative under- 
standing on H C 2 substantially. 

This paper is organized as follows. Section II consid- 
ers the weak-coupling model with gap anisotropy and 
s-wave impurity scattering. We derive an H c2 equation 
valid at all temperatures as well as an analytic expres- 
sion for H c2 (T<T c ) up to second order in 1-T/T C . The 
main analytic results of Sec. II are listed in Table I for an 
easy reference. Section III extends the H c2 equation so as 
to include p-wave impurity scattering, spin-orbit impu- 
rity scattering, and strong electron-phonon interactions. 
Section IV presents numerical examples for model Fermi 
surfaces of a three-dimensional tight-binding model with 
cubic symmetry. Section V summarizes the paper. We 
put A)b = 1 throughout. 



II. WEAK-COUPLING H c2 EQUATION 

A. Fermi-surface harmonics and gap anisotropy 

We first specify the gap anisotropy in our consideration 
with respect to the Fermi-surface harmonics. The Fermi- 
surface harmonics were introduced by Allen 58 as conve- 
nient polynomials in solving the Boltzmann and Eliash- 
berg equations. They were later used by Langmann 33 to 
derive an H C 2 equation applicable to anisotropic Fermi 
surfaces and anisotropic pairing interactions. However, 
the polynomials constructed by Allen based on the Gram- 
Schmidt orthonormalization are not very convenient for 
treating the gap anisotropy. We here adopt an alter- 
native construction starting from the pairing interaction 
y(kp,kp) on the Fermi surface, 59 where kp denotes the 
Fermi wavevector. Evidently y(kp,k F ) is Hcrmitian 
F*(kF,k F ) = F(k F ,kF), and invariant under every sym- 
metry operation R of the group G for the relevant crystal 
as RV(\tF, k F )i? _1 =y(kp,k F ). We hence consider the 
following eigenvalue problem: 

|dS F p(k F )F(k F ,k F )^)(k F ) = ^)^)(k F ). (2) 
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TABLE I: Equation numbers for the relevant analytic expressions to calculate H C 2- The upper critical field H C 2 corresponds 
to the point where the smallest eigenvalue of the Hermitian matrix A=(A N n') takes zero. 



(...) 0(k F ) 


B 


$o 


T c 




VF+ 


Cl,2 


Xij 


B 


Ann' 


K-nn' V( x ) H c2 (T<T c ) 


Bi 


B 2 




(4) (5) 


(7) 


hc/2e 


(A5) 


(15) 


(18) 


(19) 


(20) (13) 


(30) 


(36) 


(39) (41) (22) 


(23) 


(A9a) 


(24) (A4) 



Here dS F denotes an infinitesimal area on the Fermi sur- 
face and p(k F )EE [(27r) 3 Ar(0)|v F |]- 1 with v F the Fermi ve- 
locity and N (0) the density of states per one spin and per 
unit volume at the Fermi energy in the normal state. The 
superscript T denotes an irreducible representation of G, 
j distinguishes different eigenvalues belonging to T, and 
7 specifies an eigenvector in (r, j). This eigenvalue prob- 
lem was also considered by Pokrovskii 60 without specify- 
ing the symmetry. The basis functions thereby obtained 
naturally have all the properties of Fermi-surface har- 
monics introduced by Allen. Especially, they satisfy the 
orthonormality and completeness: 

(4><? i) '<t>§' i ' ) )=&rr>5 jj ,5 rr ,, (3a) 

^^wew ^y , (3b) 

where (• • • ) denotes the Fermi-surface average: 

(A) = jdS FP (k F )A(k F ). (4) 

Using Eqs. (2) and (3), we obtain an alternative expres- 
sion for the dimensionless pairing interaction A(k F , k F ) = 
-N(0)V(k F ,k' F ) as 

A(k F ,k F ) = ^A«)^)(k F )^>(k F ). (5) 

Thus, it is always possible to express a general pairing 
interaction as a sum of separable interactions. Notice 
that the above procedure is applicable also to multiband 
superconductors. Indeed, we only have to extend the 
integration over k F to all the Fermi surfaces. 

The Fermi-surface harmonics can be constructed also 
from the coupling function A(k F , k F , e n — e' n )— /i*(k F , k F ) 
in the strong-coupling Eliashberg theory, 61 ' 62 where e n = 
(2n + 1)ttT is the Matsubara energy. Indeed, we only 
have to specify an appropriate bosonic Matsubara en- 
ergy loi = 2lirT and set V"(k F ,k F ) = — [A(k F , k F ,u>i) — 
H*(k F , k' F )]/N(0) in Eqs. (2) and (3). We thereby obtain 
an alternative expression for the coupling function as 

A(k F , k F ,e„-e^)-^*(k F ,k F ) 
= £[ A ( r %„- £ 'J- M *«) ] 0f)(k F )^>(k F ) . (6) 

We expect that this construction does not depend on the 
choice of loi substantially. It is worth noting that ab ini- 
tio calculations of the coupling function are now possible 



for phonon- mediated superconductors, as performed re- 
cently for MgB 2 . 63 Hence ab initio constructions of the 
Fermi-surface harmonics by Eq. (2) can be carried out in 
principle. 

From now on we consider the cases where (i) the system 
has inversion symmetry and (ii) a single A^ rj ^ is relevant 
which belongs to an even-parity one-dimensional repre- 
sentation r. Indeed, these conditions are met for most 
superconductors. Hereafter we will drop all the indices 
as 0^ r ^(k F ) — > </>(k F ), for example, and choose <p(k F ) as 
a real function. 



B. Eilenberger equations 

Now, let us derive an H C 2 equation for the second-order 
transition in the weak-coupling model with s-wave im- 
purity scattering based on the quasiclassical Eilenberger 
equations. 64-66 The Eilenberger equations are derived 
from the Gor'kov equations by assuming a constant den- 
sity of states near the Fermi energy in the normal state 
and integrating out an irrelevant energy variable. 64-66 
Thus, phenomena closely connected with either the en- 
ergy dependence of the density of states 26 or the dis- 
creteness in the quasiparticle energy levels 13,32 ' 35-38 are 
beyond the scope of the present consideration. We also 
do not consider Josephson vortices appearing in very 
anisotropic layered superconductors. 67 Within the lim- 
itations, however, the Eilenberger equations provide one 
of the most convenient starting points for deriving an H c2 
equation, as seen below. This approach was also adopted 
by Rieck et al. 30 ' 31 

We take the external magnetic field H along the z 
axis. In the presence of Pauli paramagnetism, the av- 
erage flux density B in the bulk is connected with H as 
H = B — 4:irXnB, where Xn is the normal-state spin sus- 
ceptibility. The fact that \n i s multiplied by B rather 
than H corresponds to the fact that the spins respond to 
the true magnetic field in the bulk. It hence follows that 
B is enhanced over H as 

B = H/(l-47rx„). (7) 

The vector potential in the bulk at H = H c2 can be writ- 
ten accordingly as 

A(r) = (0,5a;, 0). (8) 

The field H is supposed to be along the direction 
(sin 9 cos ip, sin 9 sin ip, cos 9) in the crystallographic coor- 
dinates (A", Y, Z) . The two coordinate systems are con- 
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nected by the rotation matrix 



K 



cos 9 cos ip cos 9 sin ip — sin 9 

— sin f cos <p 
sin 9 cos <p sin sin ip cos 6* 



(9) 



as 1ZH = (0,0, H) T , where T denotes transpose. We as- 
sume that the vortex lattice is uniform along z. 

With the gap anisotropy specified by <f>(k F ) and in the 
presence of Pauli paramagnetism, the Eilenberger equa- 
tions read 



~-n — i/J-bB + ^- (g) + ^ hvp ■ 8 ) / 



0A+-</» 5) 



(10a) 



rp CO 

A(r)ln^ =7 rT £ 



A(r) 



(0(k F )/(e„,k F ,r)) 



(10b) 

Here /ljb is the Bohr magneton, r is the relaxation time 
by nonmagnetic impurity scattering in the second-Born 
approximation, A(r) is the pair potential, and d is de- 
fined by 



(11) 



with 4>o = hc/2e the flux quantum. We will consider 
positively charged particles following the convention; the 
case of electrons can be obtained directly by A — > —A, 
i.e., reversing the magnetic-field direction. The quasi- 
classical Green's functions / and g are connected by g = 
(l-//t)V2 sgn ( en ) with /t(e n ,k F ,r) =/*(-£„, k F ,r), 68 
and T c q denotes the transition temperature in the clean 
limit t = oo. 

To obtain B c2 , we formally expand the quasiclassical 
Green's functions up to the first order in A as / = /W 
and <7 = sgn(£„). Substituting the expressions into Eqs. 
(10a) and (10b), we obtain the linearized self-consistency 
equations as 



-/ , sgn(e„) 

e n H hvF-o 



/ (i) =0A +^-(/( 1 )), (12a) 



with 



4 = £n - i/U B -Bsgn(e„) , s n = \s n \ + — . (13) 

It 



C. Operators and basis functions 

It is useful to transform the gradient operator in Eq. 
(10a) as 



v F d — (v F+ a - v F+ (J)/V2l c ■ 



(14) 



Here l c denotes -j= times the magnetic length as 

l c = v /*o/2tt j B = yfhc/2eB . 
The operators a and are defined by 



a 


lc 


ci ic 2 




d x 


v_ 


V2 






Oy 



where the constants Ci and c 2 are constrained by 

C1C2 + c\c 2 = 2, 
so that [a, at] = l. Finally, vf+ is defined by 
v F + = c 2 v Fx + iciv Fy . 



(15) 

(16) 

(17) 
(18) 



The constants (ci,c 2 ) can be fixed conveniently by re- 
quiring that the gradient term in the Ginzburg-Landau 
equation be expressed in terms of a^a without using aa 
and ata\ i.e., the pair potential near T c be described 
in terms of the lowest Landau level only. As shown in 
Appendix A, this condition yields 



ci 



X 2 



C2 



Xy 



A/ "V — V 

\ XxxXyy Xxy 1 



XxxXyy Xxy. 
1/4 

exp itan -1 



1/4 



(19a) 



Xxy 



XxxXyy 



V 

A-xy J 



where Xij = Xioi T c) is defined by 



Xij 



24(ttT c ) s 



'w Fl u F:) -; 



(19b) 

(4>){(j)v ¥l v Fj ) 



7C(3)(« F )£lel 



(4>)(<t>VFiVFj) ((t>)(<t>)(VFiVFj) 



2re„ 



2re„ 



(2re r . 



(20) 



with £ the Riemann zeta function. Notice that \ij is di- 
mensionless, approaching to Sij as t^oo for the spher- 
ical Fermi surface. It is a direct generalization of the x 
function introduced by Gor'kov 69 to anisotropic systems. 

The operators in Eq. (16) extends (a_,a + ) introduced 
by Helfand and Werthamer 6 for anisotropic crystals. For 
uniaxial crystals, they reduce to the operators used by 
Takanaka. 18 

Using Eq. (16), we can also make up a set of basis 
functions to describe vortex-lattice structures as 70 



^Afq(r) = 



2ttZ c 



cia 2 y/irV 



JVf/2 

2 exp 

-M/2+1 



iq v \y + 



1% 



x exp 



x exp 



.nau 



1 



VWw. 



cic 2 fx - l\q y - nau 



c\l c 
]q y - nau 



(21) 
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Here N = 0, 1, 2, • • • denotes the Landau level, q is an 
arbitrary chosen magnetic Bloch vector characterizing 
the broken translational symmetry of the vortex lattice 
and specifying the core locations, and V is the volume 
of the system. The quantities a\ x and a 2 are the com- 
ponents of the basic vectors ai and & 2 in the xy plane, 
respectively, with a 2 || y and a\ x a 2 = 27rZ 2 , A/" f 2 denotes 
the number of the flux quantum in the system, and 

Hn(x) = e x2 (—■£-) e ~ x2 is the Hermite polynomial. The 
basis functions are both orthonormal and complete, sat- 
isfying a^jv q = V^V'iV-iq and a^jvq = V W+ltpN+iq- 

The function (21) is a direct generalization of the Eilcn- 
berger function 10 ^jv(r|ro) with c\ = c 2 = 1 to anisotropic 
Fermi surfaces and energy gaps. For q = in the clean 
limit, Eq. (21) reduces to the function obtained by Ricck 
et al., 30 ' 31 ' 71 However, they derived it without recource 
to the creation and annihilation operators of Eq. (16). 
These operators have simplified the derivation of the ba- 
sis functions and will also make the whole calculations 
below much easier and transparent. 



so that 

' ( v Fx) = (( v fx) cos2 <P + ( v fy) sin2 <fi) c °s 2 
+ (v FZ )sm 2 9 

(v Fy ) = (v F x) sin2 f + ( v fy) cos V 
(v Fx v Fy ) = {(v FY ) - (uf x )) cos cost/? sin 

The quantities (4>v Fx v Fy ) and (<j> v Fx VFy) can be ex- 
pressed similarly in the crystallographic coordinates once 
0(kp) is given explicitly. In particular, when 0(kp) be- 
longs to the Ai g representation, the expressions for the 
two averages are essentially the same as Eq. (26). From 
Eqs. (23), (20), and (26), we realize immediately that 
the initial slope is isotropic when (i) </>(kp) belongs to 
Ai g and (ii) the crystal has cubic symmetry. 

The expression for B 2 is more complicated as given ex- 
plicitly by Eq. (A9a). It includes Fermi-surface averages 
oiv Fx , v Fx v Fy , etc., and enables us to estimate the initial 
curvature of H c2 given the Fermi-surface structure. 



D. Analytic expression of H& near T c 



E. H C 2 equation 



Using Eq. (16), it is also possible to obtain an ana- 
lytic expression for B c2 = H c2 /(l — inxn) near T c . Let us 
express it as 



B c2 = B 1 {l-t)+B 2 (l-tf 



(22) 



with t = T/T c . The coefficients B\ and B 2 determine the 
initial slope and the curvature, respectively. 

It is shown in Appendix A that B\ is obtained as 



Bx 



247ri?$ 



7C(3)(X-X TO -X 2 y ) 1/2 (h(v F y/yT c ) 2 ' 



(23) 



where ( is the Riemann zeta function, Xij is given by Eq. 
(20), and R is defined by 



2 - 



(24) 



The factor fi(u|) 1/2 /T c in the denominator of Eq. (23) 
is essentially the BCS coherence length. 44 Also, R is di- 
mcnsionless and approaches unity for t — > oo. Equa- 
tion (23) is a direct generalization of the result by Rieck 
and Scharnberg 30 for <fi(k F ) = 1 to the cases with gap 
anisotropy and for arbitrary strength of the impurity 
scattering. 

It is convenient to express (v F iVFj) in Eq. (20) with 
respect to the crystallographic coordinates (X, Y, Z) to 
see the anisotropy in B\ manifestly. Using Eq. (9), v Fx 
and Vf v are rewritten as 



v Fx — vfx cos 9 cos ip + v F y cos 9 sin ip — vpz sin 9 
vf v — ~vfx sin ip + v F y cos ip 



(25) 



We now derive an H c2 equation which can be solved 
efficiently at all temperatures. To this end, we transform 
Eqs. (12a) and (12b) into algebraic equations by expand- 
ing A and / W in the basis functions of of Eq. (21) as 41 ' 70 



A(r) = W J2 A n VWr) • 



(27a) 



N=0 



/« (e n , k F , r) = W J2 fx (£n, k F ) VWr) • (27b) 



N=0 



Let us substitute Eqs. (14) and (27) into Eqs. (12a) and 
(12b), multiply them by tp%q(v), and perform integra- 
tions over r. Equations (12a) and (12b) are thereby 
transformed into 



^M NN ,f$=<l>* N + —(f$ ) ), 



N' 



(28a) 



A. In T -f = -,T £ (<< } >-^), (28b) 



where the matrix M is tridiagonal as 



MnN> = S' n SNN' + VN+ip*S N ,N'-i - VN/3S N ,N' + 1 , 

(29) 

with 



7w F +sgn(e n ) 
2V2l c 



(30) 
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We first focus on Eq. (28a) and introduce the matrix 
/Cby 

K NNI = (■M~ 1 )nn i , (31) 
which necessarily has the same symmetry as M: 72 



K-NN'i^m P) — /CjV'Jv( e rai — /?*) — K-NN'(~ £ n, ~ P*) 



Using JC, Eq. (28a) is solved formally as 

fV = -£ic NN , (Vw + A</«) 



AT' 



(32) 



(33) 



Taking the Fermi-surface average to obtain (fffl) and 
substituting it back into Eq. (33), we arrive at an ex- 
pression for the vector fW = (f^, fi~\ • • • ) T as 



f (i) 



K.(j)+—!C [1 
2r 



h(ic) 

2r 



(34) 



with I the unit matrix in the Landau-level indices and 

A^(A ,A 1 ,A 2 ,-..) T . 

We next substitute Eq. (34) into Eq. (28b). We thereby 
obtain the condition that Eq. (28b) has a nontrivial so- 
lution for A as 



dct A = , 

where the matrix A is defined by 



(35) 



A = l\n— + ttT J2 



l cO 



(JC<f>) 



(36) 



with X the unit matrix in the Landau-level indices. The 
upper critical field B c2 corresponds to the highest field 
where Eq. (35) is satisfied, with B and H connected by 
Eq. (7). Put it another way, B c2 is determined by re- 
quiring that the smallest eigenvalue of A be zero. Notice 
that A is Hermitian, as can be shown by using Eq. (32), 
so that it can be diagonalized easily. 

Equation (36) tells us that central to determining B c2 
lies the calculation of JCnn' defined by Eqs. (29) and 
(31). An efficient algorithm for it was already developed 
in Sec. IIF of Rcf. 41, which is summarized as follows. 
Let us define R N {N = 0, 1,2, • • •) and R N (AT = 1,2, • • •) 

by 



fljv-i = (1 + Nx 2 R N )~ 1 



R N+1 = (I + NxZRn)- 1 
respectively, with 

x=\P\/? n 



Ri = l, 



(37a) 



(37b) 



(38) 



Then JCmn' for N>N' can be obtained by 

N-N' 



/Cjvjv = -^-rjN{x)riN>{x) ( — 



with 



N 



m = ViV! Y[ R k , 



fe=0 



Vn = < 



(N = 0) 



N 



(39) 



(40a) 



(40b) 



The expression of JC nn> for N <N' follows immediately 
by Eq. (32). 

As shown in Appendix B, Eqs. (40a) and (40b) can be 
written alternatively as 



2 r°° s N 



1 



s n _Hn(s) 
+ 2x 2 s 2 



ds 



'Nl 



s N exp ( — s — ^rs 2 ) ds 



2 N V^W.z N+1 e z2 i N erk(z), 



(41a) 



7] N (x) 



N 



H 



N 



Nl \V2iJ 
yVM \dy. 



V2x 

N 



y 2 /2 



(41b) 



y=l/x 



respectively, where z=l/V2x and i Jv erfc(z) denotes the 
repeated integral of the error function. 73 The latter func- 
tion t}n(x) is an ^th-order (^3-th-order) polynomial of 
x 2 for A^ = cvcn (odd). 

Thus, the key quantity /CjviV' is given here in a com- 
pact separable form with respect to N and N'. This is a 
plausible feature for performing numerical calculations, 
which may be considered as one of the main advantages 
of the present formalism over that of Langmann. 33 Our 
/Coo in Eq. (39) is more convenient than Eq. (26) of Ho- 
henberg and Werthamer 9 in that H c2 near T c is described 
in terms of the lowest Landau level for arbitrary crystal 
structures. 

Equations (35) and (36) with Eqs. (39), (13), (30), 
(15), (18), (20), and (19) are one of the main results of 
the paper (see also Table I). They enable us efficient cal- 
culations of H c2 at all temperatures based on the Fermi 
surfaces from ab initio electronic-structure calculations. 
They form a direct extension of the Rieck-Scharnbcrg 
equation 30 to the cases with gap anisotropy and arbi- 
trary strength of the impurity scattering. Indeed, Eq. 
(41b) is written alternatively as 



V2N(x) 



1 



y/(2Wy.2 N z 



2..Y 



(42) 
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with z = l/\/2x, where Pn is the polynomial defined 
below Eq. (6) of Rieck and Scharnberg. 30 Substituting 
this result and the last expression of Eq. (41a) into Eq. 
(39), it can be checked directly that e n K,2N'2N for N' <N 
is equal to M 2 n>2N in Eq. (6) of Rieck and Scharnberg. 30 
Using this fact, one can show that the matrix A in Eq. 
(36) reduces to the corresponding matrix in Eq. (5) of 
Rieck and Scharnberg either (i) for the isotropic gap with 
arbitrary impurity scattering or (ii) in the clean limit 
with an arbitrary gap structure. Here we have adopted x 
in Eq. (38) as a variable instead of z, because x remains 
finite at finite temperatures. 

From Eq. (39) and the symmetry (3 — > — /3 for vp — > 
— vf, we realize that (K,2N,2N>+i) , {fc2N,2N>+i<fr), and 
(K-2N \2N '' +i4> 2 ) all vanish in the present case where (i) 
the system has inversion symmetry and (ii) </>(kp) be- 
longs to an even-parity representation. It hence follows 
that we only have to consider TV = even Landau levels 
in the calculation of Eq. (36). To obtain a matrix ele- 
ment of Eq. (36), we have to perform a Fermi surface 
integral for each n and perform the summation over n, 
which is well within the capacity of modern computers, 
however. Actual calculations of the smallest eigenvalue 
may be performed by taking only N < N cut Landau lev- 
els into account, and the convergence can be checked by 
increasing iV cut . We can put N cut = near T c due to 
Eq. (19), and have to increase N cut as the temperature is 
lowered. However, excellent convergence is expected at 
all temperatures by choosing A cut <20. 



First of all, we derive expressions for T c at H = 0, the 
coefficients (ci, C2) in Eq. (16), and B C 2 near T c up to the 
first order in 1 — t, based on Eq. (44) and following the 
procedure in Sec. A. It turns out that we only need a 
change of the definition of Xij from Eq. (20) into 



\3 00 



Xij 



<?3 



24(7rT e 

7 C(3)(4)^ ^ 
+ 



V Fi V F j 

3 



■ + 



(0) 



2re n /h 



2T 1 e n /h 

where the matrices V and Q are defined by 

(0) 



Tij — 



<j>- 



2T\e n \/H 



3ft 



FiVFj 



2Ti£„ 



(45) 



(46a) 



(46b) 



Then T c , (01,02), and B\ in Eq. (22) are given by the 
same equations, i.e., Eqs. (A5), (19), and (23), respec- 
tively. 

Using Eqs. (14) and (27), we next transform Eq. (44) 
into an algebraic equation. The resulting equation can 
solved in the same way as Eq. (33) to yield 

f (1) = JC(*A + A (f (i)) + |Lk F . ( k F f(i))') , (47) 



III. EXTENSIONS OF THE H c2 EQUATION 

We extend the H c2 equation of Sec. II in several direc- 
tions. 



where K, is given by Eq. (39). It is convenient to introduce 
the quantities: 



Pj 



^-k Fj (j = x,y,z). (48) 



A. p-wave impurity scattering 

We first take p-wave impurity scattering into account. 
In this case, Eq. (10a) is replaced by 



ft 3/i 1 

e n -i m B + — (g) + —k F -(k / F gy + -hv F -d)f 



1 (^+| I (/) + ^-kF-(kF/) / ).'/- 



(43) 



where (k' F g)' = (k. F g(e n , k F , r))', for example, and kp = 
kp / (fcp) x l 2 . Notice that kp is not a unit vector in general. 
Linearizing Eq. (43) with respect to A, we obtain 



4 + ^^p.a)/(^0A + A (/ ( 1 ) ) 

3h. 



+ ^k F .(k F/ «)', (44) 



with e' n defined by Eq. (13). 



Then from Eq. (47), we obtain self-consistent equations 
for (p* f ( V) and (pjf* 1 )) as 



(49) 



- (p* fw) - 






(p*J w ) 
(p;f (1) > 


= w 


(p*/C0)A 


. (p*j {1) ) . 







where the matrix W is defined by 

'Z-dPoffc) -(PoPxlC) -(poPyK.) 
-(plPoK.) l-{\p x \ 2 K) -(p* x p y JC) 

-{P* V P0K) -(p* y P X K) l-(\Py\ 2 lC) 



We 



■(PoPzty ' 
■{PlPzlQ 

-{PyPzK) 

»/C) -lp* zPx K) -<p* zPy )C) l-(\ Pz \ 2 IC) 



(50) 

The complex conjugations * in Eqs. (49) and (50) are 
not necessary here but for a later convenience. Notice 
the symmetry W* m (s n , 0) = W m j (-£•„, 0) in the matrix 
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elements of W, as seen from Eq. (32). Using Eq. (49) in 
Eq. (47), we obtain an explicit expression for as 



f (1) = K,(pA + [p K p x K PylC Pz K\W 



(p* x K,4>) 



A. 



(51) 

Finally, let us substitute Eq. (51) into Eq. (28b). We 
thereby find that Eq. (36) is replaced by 



T 

A=lln — 

J-cO 



7TT 



oo 

E 



ITT - W> 



-[<PoAC0) (p x /C0) faK® (p z W)]W 



(PolC<j>) 
(p* x K,(p) 



(52) 



As before, H C 2 is determined by requiring that the small- 
est eigenvalue of Eq. (52) be zero. This A is Hermitian, 
as can be shown by using Eq. (32) and W ; * m (£„,/?) = 
W TO ;(— £„, /?). Thus, Eq. (52) can be diagonalized easily. 

It is straightforward to extend Eq. (52) to a more gen- 
eral impurity scattering with the k F -dependent relax- 
ation time r(kp,kp). To this end, we apply the pro- 
cedure of Eqs. (2)-(5) to l/r(kp,kp) to expand it as 



f 



r(k F ,k F ) 



E 

rj7 



(53) 



where l/r' r ^ and ^^(kp) denote an eigenvalue and its 
eigenfunction, respectively. We then realize that 



2r«) 

substitutes for p and pj in Eq. (52) 



(k F ) 



(54) 



B. Spin-orbit impurity scattering 

It was noticed by Werthamcr et al. 7 and Maki 8 that, 
for high-field superconducting alloys with short mean free 
paths, Pauli paramagnetism has to be incorporated si- 
multaneously with spin-orbit impurity scattering. They 
presented a theory valid for r -C t so , where r so is spin- 
orbit scattering time. It was later generalized by Rieck 
et al. 31 for an arbitrary value of r so . This effect can also 
be taken into account easily in the formulation. 

In the presence of spin-orbit impurity scattering, Eq. 
(10a) is replaced by 

e n -i m B+^(g) + ^(\k F xk F \ 2 g y + ±hv F -d\f 



U+A (/) + g^ ( |k Fxk ^|2 /r ) (/ . 



(55) 



with c so = l/((|kFxk F | 2 )'). To simplify the notations and 
make the argument transparent, it is useful to introduce 
the quantities: 



Po= V|l. PiJ = J ^ (^f<% - k F ik F j) , (56a) 



qo " V 27 ' qi ^ = \j ^ k Fi k Fj (2 - % ) , (56b) 
and the vectors: 

P = {PQ,Pxx,PyyTPzz,Pxy,Pyz,Pzx) 7 (57a) 
Qxx 1 Qyy 7 Qzz ; Qxy -> Qyz i Qzx )■ (57b) 

Then Eq. (55) linearized with respect to A is written in 
terms of Eq. (57) as 

(4 + ^W©) /« = ^A + p • (q/«) , (58) 

where e' n is defined by 

i' n = i n - i/j, B Bsgn(s n ) , i n = \e n \ + p • (q) . (59) 

Notice p • (q) = (p) • q. 

It follows from the procedure in Sec. A that T c at H = 
satisfies 



±c n=0 



£. 



where the matrix Q is defined by (r, s = 

fq r p. 



Q 



rs u rs 



(60) 

c) 

(61) 



Also, Xij m Eq- (20) should be modified into 



24(7rT e ) 3 " 
X " := 7C(3)(,|)^ o 



V Fi V Fj 



X 



■+p T e- 1 (P 



(62) 



Finally, i? in Eq. (24) is replaced by 

00 r ' A2. 



i? = 1 - 2ttT c ^ 

/ p-(q)p T ^ 



+ 



2 p-(q> 



Q- 1 ^ 



(63) 



With the above modifications, T c , (c 1? C2), and .Bi in Eq. 
(22) are given by Eqs. (A5), (19), and (23), respectively. 



9 



We now transform Eq. (58) into an algebraic equation 
by using Eqs. (14) and (27). The resulting equation can 
solved in the same way as Eq. (33). We thereby obtain 



(64) 



where K. is given by Eq. (39) with e' n replaced by Eq. 
(59). From Eq. (64), we obtain self-consistent equations 
for (q i {1) ) and (tfyf* 1 )) as 



" (5of«> " 




- (g o /C0)A " 


(fef (1) ) 




(q xx K,(j>)A 


<<^f (1) > 




(q vy )C<f>)A 


(fef (1) ) 


= W 


{q zz JC4>)A 


<^f (1) > 




(q xy JC(t>)A 


<^f (1) > 




(q yz K,<p)A 


. (fef (1) > . 




. {q zx K,(j))A _ 



(65) 



thereby find that Eq. (36) is replaced by 



w 



-[(p Q JC4>) (p xx lC4>) {p yy JC4>) 

n— — oo 

-[(qoJC<j>) (q xx JC<j>) (q yy JC(j)) ■ ■ -]W^ 



(qo£<t>) 

(q xx K,<j)} 
{q vv K,(t>) 



(p xx }C<f>) 
{p yy K,4>) 



(68) 



where the matrix W is defined by 



As before, H C 2 is determined by requiring that the small- 
est eigenvalue of Eq. (68) be zero. This A is Hcrmitian, 
as can be shown by using Eq. (32) and [W^(e n , j3)]i m = 
W^j (—£„,/?), which can be diagonalized easily. 



W 



1-{qoPofc) -(qoPxxIC) -(q p yy IC) 

-(q-xxPolC) Z-(q X xPxx)C) -(qxxPyylC) 



~(q VV PolC} -{qyyPxxIC} 

-(q zz p !C) -{qzzPxxk) 

-(qxyPolC) -{q xy p xx K) 

-WyzPofc) -{qyzPxxlC} 

-(qzxPolC) -{qzxPxxK.) 



T- (lyyPyyfc) 
-(qzzPyyIC) 

~ (ixyPyyfc) 
-(qyzPyyJC) 
-(qzxPyylC) 



(66) 

Using Eq. (65) in Eq. (64), we obtain an explicit expres- 
sion for f « as 



f (1) = K(j)A+ [poK. PxxIC PyyK 



W 



K,(j)A+ [q K q xx K, q yy K, ■ ■ ■ ] W f 



{qoK<t>) ' 

(q X xlC(j>) 

(q yy K.<t>) 

(PxxK-cj)) 

{p yy K,4>) 



(67) 



with Wt defined by [WHe n J)} lm = W* ml (-e n J). The 
latter expression originates from the self-consistency 
equations for (pof' 1 ') and (p^fW) similar to Eq. (65). 
Finally, let us substitute Eq. (67) into Eq. (28b). We 



C. Strong electron-phonon interactions 

We finally consider the effects of strong electron- 
phonon interactions within the framework of the Eliash- 
berg theory. 61 ' 62 We adopt the notations used by Allen 
and B. Mitrovic 62 except the replacement ZA^A. 

The Eilenberger equations were extended by Teichler 74 
to include the strong-coupling effects. They can also be 
derived directly from the equations given by Allen and 
B. Mitrovic 62 by carrying out the "£ integration" 66 as 

Ze n -i m B+^-{g) + \tw ¥ -ay = (A0+|-(/)) ff , 

(69a) 

A(e„,r)=7rT {\(e n -e nl )-n*](<j){k F )f(e n/ , k F , r)) , 

n' =— n c o 

(69b) 

Z(£„,k F ) = H ^ (A(k F ,k F , £„-£„,)#(£„/, k F ,r))', 

n'— — n c (i 

(69c) 

where n c o corresponds to the Matsubara frequency about 
five times as large as the Debye frequency. 62 We have 
retained full kp dependence of A in Eq. (69c), because 
the contribution from other pairing channels, which may 
be negligible for the pair potential, can be substantial for 
the renormalization factor Z. 

We linearize Eqs. (69) with respect to A and repeat the 
procedure in Sec. A up to the zeroth order in 1—t. It then 
follows that T c at H — is determined by the condition 
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that the smallest eigenvalue of the following matrix be 
zero: 



A (0) =6 , 



7rT[A(£„-e n 0-M*] 
2r\e tl 



(70) 



where is given by 



Z< >(e„,k F ) = l- 



ttT 



^2 (A(k F ,k / F ,e n -e n /)> / sgn(e n /), 

=-n c0 

(71) 



and e„ is defined together with e' n by 



e n = Z<-% n \ + £.. 

It 



i/i B -Bsgn(e„) . (72) 



We next fix (ci,c 2 ) in Eq. (16) conveniently. For the 
weak-coupling model, we have fixed it by using Eq. (A6) 
near T c so that the coefficient of aa vanishes, i.e., there 
is no mixing of higher Landau levels in the H c2 equa- 
tion near T c . However, the coefficient of aa in the cor- 
responding strong-coupling equation becomes frequency 
dependent. It hence follows that, even near T c , there 
is no choice for (ci,C2) which prevents mixing of higher 
Landau levels from the H c2 equation. We here adopt the 
weak-coupling expression in Eq. (19). 

We now consider the H c2 equation and repeat the same 
calculations as those in Sec. IIB. We thereby find that Eq. 
(36) is replaced by 



A-nN,n' 



N> 



Snn'^NN' ~ [A (e„ -£„')- £t* 



(K'tff 



(/C'</>) 



(73) 



NN' 



where KJ = K.(js n /,/3) which also has kp dependence 
through = Z(°)(e n /,k P ). As before, H c2 is deter- 
mined by requiring that the smallest eigenvalue of Eq. 
(73) be zero. 

We may alternatively use, instead of Eq. (73), the ma- 
trix: 



nN,n'N' ~ A**)nn'^W"JV 

h 



It 



(K4>) 



(Kct>) 



(74) 



NN' 



where (A — /i*) _1 denotes inverse matrix of A — fj,*. It is 
Hermitian for n^B— >0, and also acquire the property by 
combining n>0 and n<0 elements. 



IV. MODEL CALCULATIONS 

We now present results of a model calculation based on 
the formalism developed above. We restrict ourselves to 




FIG. 1: Fermi Surfaces of the tight-binding model in the sim- 
ple cubic lattice. The Fermi energies are: (a) £f = —3, (b) 
-2, (c) -1, and (d) 0. 



the weak-coupling model of Sec. II with an isotropic gap, 
no impurities, and no Pauli paramagnetism. As for the 
energy-band structure, we adopt a tight-binding model 
in the simple cubic lattice whose dispersion is given by 



£k = ~2i {cos(k x a) + cos(k y a) + cos(k z a)} . 



(75) 



Here a denotes lattice spacing of the cubic unit cell and t 
is the nearest- neighbor transfer integral. We set t = a = 1 
in the following. The corresponding Fermi surfaces are 
plotted in Fig. 1 for various values of the Fermi energy 
£p- For ep ~ —6, i.e., near the bottom of the band, the 
Fermi surface is almost spherical with slight distortion 
due to the cubic symmetry. As £f increases, the cu- 
bic distortion is gradually enhanced. Then at £f = — 2, 
the Fermi surface touches the Brillouin-zone boundary at 
k x = (0, 0, ±tt), (0, ±, 7r, 0), (±tt , 0, 0). Above this critical 
Fermi energy, the topology of the Fermi surface changes 
as shown in Fig. 1(c). It is interesting to see how such a 
topological change of the Fermi surface affects H c2 . 

We computed H c2 based on Eq. (35) in the clean limit 
without Pauli paramagnetism. The Fermi-surface aver- 
age in Eq. (36) was performed by two different methods. 
For general values of sp, we used the linear tetrahedron 
method which is applicable to any structure of the Fermi 
surface. In this method, the irreducible Brillouin zone 
is divided into a collection of small tetrahedra. From 
each tetrahedron which intersects the Fermi surface, a 
segment of the Fermi surface is obtained as a polygon 
by a linear interpolation of the energy band. Numerical 
integrations over the Fermi surface were then performed 
as a sum over those polygons. Another description of the 
Fermi surface is possible for £p < — 2, where we can adopt 
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TABLE II: The ratio B 2 /Bi for the field directions [100], 
[110], and [111] in the cases £f = —3 and —2.02. The quantities 
B\ and B2 are defined in Eq. (22). The values should be 
compared with 0.13 for the spherical Fermi surface. 



0.0 0.2 0.4 0.6 0.8 1.0 

T/T r 



FIG. 2: Curves of the reduced critical field h* d (i) for the cubic 
tight-binding model with £f = —2.02 (dotted lines), £f = —3 
(solid lines), and e F — > —6 (i.e., the spherical Fermi surface; 
dash-dotted line). The field directions are d= [111], [HO], and 
[100] from top to bottom in each case. 



the polar coordinate k = (k sin 9 cos 0, k sin 9 sin </>, k cos 9) 
and the Fermi surface kp = kp(9, <p) is obtained by solv- 
ing the equation £k = £f numerically for each (9,(j)). 
An integration over the Fermi surface is then performed 
by using the variables (9,(j)). We performed both types 
of calculations to check the numerical convergence of the 
tetrahedron method. Excellent agreements were achieved 
generally by using 3375 tetrahedrons. An exception is the 
region £p « —2, where larger number of tetrahedrons was 
necessary due to the singularity around kx- 

The infinite matrix Ann' in Eq. (36) was approxi- 
mated by a finite matrix of N, N' < N cut , and the con- 
vergence was checked by increasing N cut . The choice 
-^cut = is sufficient for T<T C , and it was found numeri- 
cally that N cut = 8 yields enough convergence for all field 
directions at the lowest temperatures. It was also found 
that higher Landau levels of N > 1 contribute to H c2 by 
only 4% even at T/T c = 0.05. Thus, the lowest-Landau- 
level approximation to the pair potential is excellent for 
this cubic lattice. This is not generally the case, however, 
and the contribution of higher Landau levels can be con- 
siderable for low-symmetry crystals, as will be reported 
elsewhere. 75 

Before presenting any detailed results, it is worth not- 
ing that the GL equations, 1,2 where the anisotropy en- 
ters only through the effective-mass tensor, cannot ex- 
plain possible anisotropy of H c2 in cubic symmetry, as 
already pointed out by Hohcnberg and Werthamer. 9 This 
GL theory is valid near T c so that the upper critical field 
for T<T C should be isotropic in the present model. The 
anisotropy of H C 2 in cubic symmetry emerges gradually 
at lower temperatures, as seen below. 



EF 


[100] 


[110] 


[111] 


-3 


0.08 


0.27 


0.33 


-2.02 


0.44 


0.78 


0.90 



We calculated the reduced critical field h* (t) defined by 
Eq. (1) for the magnetic field directions d = [100], [110], 
and [111]; we denote them as h* d (t). Figure 2 presents 
h* d {t) for £ F = -3 and -2.02 as a function of t = T/T c . For 
£f = —3, h*(t) is almost isotropic for t > 0.8 and cannot 
be distinguished from the curve for the spherical Fermi 
surface. At lower temperatures, the anisotropy appears 
gradually. Whereas hf 1QO i (t) is reduced from the value 
for the spherical Fermi surface, /i* m ](i) and h* 11Q ^(t) are 
enhanced due to the cubic distortion of the Fermi sur- 
face. At t — 0.05, ^* m ](i) and h* lw ^(t) are larger than 
/i* 100 j(i) by 19% and 15%, respectively. In another case 
£f = —2.02 where the Fermi surface nearly touches the 
Brillouin zone boundary, h* d {t) are remarkably enhanced 
for all field directions. Especially, h* ul ^(t) and ft* 110 ](i) 
at low temperatures exhibit values about 60-70% larger 
than those for the spherical Fermi surface. 

At £p = —3, /i* 111 ](i) and h^ 11Q ^(t) near T c show small 
upward curvature, whereas ht QQ , (t) remains almost iden- 
tical with the curve for the spherical Fermi surface. This 
difference may be quantified by the ratio B 2 /Bi defined 
in Eq. (22). It was numerically evaluated by using the 
Fermi velocity on the Fermi surface and shown in Ta- 
ble II. The values for the directions [110] and [111] are 
larger than 0.13 for the spherical Fermi surface. Thus, 
calculated B 2 /Bi values well describe the difference in 
h* (t) for t < 1 among field directions. The upward cur- 
vature is more and more pronounced as the Fermi sur- 
face approaches the Brillouin zone boundary, as can be 
seen clearly in Fig. 2 for £p = —2.02. The correspond- 
ing ratio B2/B1 for the [110] and [111] directions are 
about three times larger than those for £p = —3. Thus, 
the present calculation clearly indicates that the Fermi 
surface anisotropy can be a main source of the upward 
curvature in H C 2 near T c . 

In Fig. 3, we plot h* d {t) at £ = 0.05 as a function of £ F . 
As £p^— 6, the angle dependence of h* d (t) vanishes and 
it converges to the value for the spherical Fermi surface. 
As £p is increased from —6, cubic distortion is gradually 
introduced to the Fermi surface as shown in Fig. 1, and 
h* d {t) gradually develops anisotropy as a consequence. 
For — 6<£f<— 2.5, curves of h* 100 ^(t) fall below that for 
the spherical Fermi surface, whereas /im 10 ] (t) and (t) 
are enhanced over it. As £p approaches to —2, h* d {t) is 
enhanced significantly irrespective of the field direction. 
Indeed, h* d (t) for every field direction shows a singularity 
at £f = — 2 where the Fermi surface touches the Brillouin 
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0.10 



FIG. 3: The reduced upper critical field h* d {t) at t = 0.05 as 
a function of the Fermi energy sf- The field directions are 
d= [111], [110], and [100] from top to bottom, respectively. 



zone at kx with vanishing Fermi velocity vp at these 
points. As a result, the contribution around these points 
becomes important in the integration (Knni) over the 
Fermi surface at low temperatures. This is the origin of 
the enhancement of h* d (t) around £f = — 2. For e-p > — 2, 
the difference between /i* 110 j and is larger than that 
for ep<— 2.5. This may be attributed to the topological 
difference of the Fermi surface. At £p = 0, the tight- 
binding band is half-filled and the Fermi-surface nesting 
occurs. However, h* d (t) docs not show any singularity 
around this energy. 

Finally, we present results on the higher Landau-level 
contributions to the pair potential A(r) which is ex- 
panded as Eq. (27a). In general, when the system has 
n-fold symmetry around the field direction, mixing of 
higher Landau levels with multiples of n develops as 
the temperature is lowered. 70 Figure 4 shows the ratio 
Ajv/Aq as a function of T/T c for ep = —3 (solid lines) 
and £ F = -2.02 (dotted lines) with (a) H|| [100] {N = 4,8 
from bottom to top lines), (b) H|| [110] (iV = 2,4,6 from 
bottom to top lines), and (c) H | [111] (N = 6). One 
can clearly observe a general tendency that the mixing 
is more pronounced as the symmetry around H becomes 
lower as well as £p approaches closer to —2. Especially 
when H || [110] and e F = -2.02, the N = 2 contribution 
reaches up to nearly 15% of the lowest Landau-level con- 
tribution as T^0. The results suggest that the lowest- 
Landau-level approximation for the pair potential 9 is not 
quantitatively reliable at low temperatures for the field 
along low-symmetry directions, for complicated Fermi 
surfaces with divergences in the components of vp per- 
pendicular to H, or for low-symmetry crystals. 



V. SUMMARY 

We have derived an efficient H c2 equation incorporat- 
ing Fermi-surface anisotropy, gap anisotropy, and impu- 
rity scattering simultaneously. Basic results of Sec. II are 
summarized in Table I. This H C 2 equation is a direct ex- 
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FIG. 4: The ratio Ajv/Ao of the expansion coefficients in 
Eq. (27a) as a function of temperature with (a) H | [100], 
(b) H || [110], and (c) H || [111]. The solid and dotted lines 
correspond to £f = —3 and £f = —2.02, respectively, with (a) 
N = 4, 8 from bottom to top (b) N = 2, 4, 6 from bottom to 
top, and (c) iV = 6. 



tension of the Ricck-Scharnberg equation 30 and reduces 
to the latter cither (i) for the isotropic gap with arbi- 
trary impurity scattering or (ii) in the clean limit with 
an arbitrary gap structure, as shown around Eq. (42). 
The operators introduced in Eq. (16) have been help- 
ful to make the derivation simpler than that by Rieck 
et al. 3a ' 31 The present method will be more suitable for 
extending the consideration to multi-component-order- 
parameter systems or to fields below H c2 . 

We have also obtained a couple of analytic expressions 
near T c (i) for H c2 up to the second order in 1 — T/T c 
and (ii) for the pair potential up to the first order in 
1-T/T C . The latter result is given by Eq. (A8) with Eqs. 
(A9b) and (A4). They are useful to estimate the initial 
curvature of H C 2 as well as the mixing of higher Landau 
levels in the pair potential. 
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The H C 2 equation of Sec. II has also been extended in 
Sec. Ill to include p-wave impurity scattering, spin-orbit 
impurity scattering, and strong electron-phonon interac- 
tions. 

Finally, we have presented numerical examples in Sec. 
IV performed for model Fermi surfaces from the three- 
dimensional tight-binding model. The results clearly 
demonstrate crucial importance of including detailed 
Fermi-surface structures in the calculation of H c2 . It 
has been found that, as the Fermi surface approaches 
the Brillouin zone boundary, the reduced critical field 
h*(t) in Eq. (1) is much enhanced over the value for the 
isotropic model with a significant upward curvature near 
T 

It is very interesting to see to what degree the up- 
per critical field of classic type-II superconductors can 
be described quantitatively by calculations using realistic 
Fermi surfaces. The result by Butler 22,23 on high-purity 
Niobium provides promise to this issue. We have per- 
formed detailed evaluations of H C 2 for various materials 
based on Eq. (35) by using Fermi surfaces from density- 
functional electronic-structure calculations as an input. 
The results are reported elsewhere. 75 
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APPENDIX A: DETERMINATION OF (ci,c 2 ) 
AND ANALYTIC EXPRESSION OF H c2 NEAR T c 

We here fix the constants (ci, c 2 ) in Eqs. (16)-(18) con- 
veniently so that H C 2 near T c can be described in terms 
of the lowest Landau level only. We also derive analytic 
expressions for B\ and B 2 in Eq. (22) so that one can 
calculate them once the relevant Fermi-surface structure 
is given. 

In the region T <T C where l c ^>oo in Eq. (14), we can 
perform a perturbation expansion with respect to the 
gradient operator v F c?. The equation for the i^th-ordcr 
solution f^p (f = 0, 1, • • • ) is obtained from Eq. (12a) as 



4>A h{fl 1] ) sgn(e„) (1) 



with /IV =0. Noting 0(-k F ) = 0(k F ), we solve Eq. (Al) 
self-consistently for (fi^), put the resulting expression 
back into Eq. (Al) to express explicitly, and finally 
take the Fermi-surface average {(pf'tP)- This procedure 
yields 



m 2 

2r\s n \ 



(A2a) 



W 2 (1) ) - 



Wi 1} > 





(av F • d) 2 ) A , (A2b) 



with \e n \' = \s n \-iiJ,BBsga(e n ), and ((f)/^) = ((pf^) =0. 

Let us substitute Eq. (A2) into Eq. (12b), replace the 
gradient operator by the right-hand side of Eq. (14), put 
B = B c2 in l c of Eq. (15), and expand \e n \' with respect 
to HBB C 2/\e n \. We thereby obtain the self-consistency 
equation near T c as 

WQfiA + \_W 2 .2 O^o) + W2.2 aa ~ W2fii aa ^ +0* a )] A 

+ [ — — J [w 4 4 a) a) a) a) + w\ 4 aaaa 
V Bi J 

—w 4,2(01' a) a' + a) aa) a) + a) a) aa) + a) a) a) a) 
—w\ 2 (a' aaa + aa) aa + aaa) a + aaaa^) 
+!U4.0o(«a aa' + a' aaa' + aa)a) a + a^aa'a) 
+Wifib(aaa) a) + a^a^aa) + top] A = . (A3) 

Here B\ is given in Eq. (22), which is incorporated into 
the denominator for convenience. The functions w v ,n = 
w u p (T) and wp = wp (T) are dimensionless and defined 

by' 

W0fi (T) ee In ^ - (1- {^) 2nT £ ( 1 - 1 



n=0 



(A4a) 



W 2 .2(T) 



W2,o(T) 



BihWT 




, (A4d) 



BjhWT ^ 1 

n=0 ™ 



8$ 2 
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2re„ 




2reJ Wf+|Wf+I 



2re„/ 



VF+\ 



(A4e) 
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w P {T) = -(mbSx) 2 2ttT^ 
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Q "I" 



(A4g) 
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We next substitute Eq. (22) into Eq. (A3) and expand 
w„ tll in Eq. (A3) up to the ^^th order in 1 — t. We 
also put wp(T) — wp(T c ). This procedure yields three 
equations corresponding to order 1, 1 — t, and (1 — t) 2 . 
The equation of order 1 is given by wo,o(7c)A = 0. It 
determines T c at H — Q by 



In 



T, 



cO 



l h 

' ; ' 2 47rrT f 



(A5) 



with ip(x) the digamma function. 

The equation of order 1 — t in Eq. (A3) is obtained as 

[-T C < (T C ) - w 2fi (T c )(2tJa+l) 

+u>2,2(T c )a t a t + w^ 2 (T c )aa] A(r) = . (A6) 

To solve it, we use the arbitrariness in (ci, c 2 ) and impose 
W2,2(T C ) = 0. Noting Eqs. (A4b) and (18), this condition 
is transformed into a dimensionlcss form as 



XxxC 2 + 2iXxyCiC 2 



Xyy c l 



0. 



(A7) 



where Xij = Xij(T c ) is defined by Eq. (20). Equation 
(A7) can be solved easily in terms of c 2 . Substituting 
the resulting expression into Eq. (17) and choosing c\ 
real, we obtain Eq. (19). 

Now that w 2 , 2 (T c ) = in Eq. (A6), the highest field for 
a nontrivial solution corresponds to the lowest Landau 
level where w 2 fi(T c ) = —T c w' 00 (T c ). Introducing R = 
—T c w' (T c ) which is given explicitly as Eq. (24), and 
using Eqs. (A4a), (A4c), (18), and (19), we obtain the 
expression for B\ as Eq. (23). 

We finally consider the equation of order (1— t) 2 in Eq. 
(A3) and expand the pair potential as 

A(r) = A {Vo q (r) + (l-t)[r 2 ^ 2q (r)+r 4 ^4 q (r)]} , (A8) 

where ?/Vq( r ) is defined by Eq. (21), and (A ,r 2 ,r 4 ) are 
the expansion coefficients with (r 2 , r 4 ) describing relative 
mixing of higher Landau levels in the pair potential. Let 



us substitute Eq. (A8) into Eq. (A3), multiply the equa- 
tion of order (1 — t) 2 by ^^ q ( r )j an d perform integration 
over r. The resulting equations for N = 0, 2, 4 yield 



B, = 



'l+T c w' 20 +w i fi a +2w4 t ob+wp 



R 



T C W 2 2 + 6w 4 ,2 

r 2 = '—F= 

2y/2R 



\Z6u>4 4 

ri = ^R- 



Bi, (A9a) 



(A9b) 



respectively. The functions in Eqs. (A9a) and (A9b) are 
defined by Eqs. (A4) and (24) and should be evaluated 
at T c . In the clean limit r^oo, these functions acquire 
simple expressions as 

R = T?w$ ) = l, 7>'=-2, T c w'=0, (AlOa) 



W4,n = 



31C(5) 
[7C(3)P 



\VF- 



|4-M, 



2\2 



Wp 



(ct> 2 \v F+ \ 2 ) 

7C(3)( MB i?i) 2 
4(^T C ) 2 



(AlOb) 



(AlOc) 



with \x = 0,2,4 and 11)4,0 = u> 4 ,oa — Wi.ob- Equation 
(A9) with Eq. (A10) includes the result by Hohenberg 
and Wcrthamcr 9 for cubic materials, and also the one 
by Takanaka 18 for uniaxial materials in the relevant or- 
der, both except the Pauli term wp. Thus, we have ex- 
tended the results by Hohenberg and Werthamer 9 and by 
Takanaka 18 to arbitrary crystal structures and impurity- 
scattering time, including also Pauli paramagnetism. 

Equation (A9) reveals a close connection of both the 
curvature in H c2 (T <T C ) and the mixing of higher Lan- 
dau levels in A(r) with the Fermi-surface structure. For 
example, we realize from Eq. (A9b) with Eqs. (A10), 
(18), and (26) that the mixing of N = 2 Landau level 
is absent for cubic materials where c\ = c 2 = 1. This is 
not the case for low symmetry crystals, however. Equa- 
tion (A9) enables us to estimate the curvature and the 
mixing based on Fermi-surface structures from detailed 
electronic-structure calculations . 



APPENDIX B: PROOF OF EQ. (41) 

The first expression in Eq. (41a) can be proved by in- 
duction as follows. First of all, r] = R is transformed 
from Eq. (37a) as 76 



Vo 



V2 



l/2x' 



1 + 



/•OO 

Jl/s/2x 



c s ds 



1 + 



2x z 



1 + 



f cxp {~ s -i s2 ) ds = 7*L T 



+ 2x 2 s 2 



ds 



00 2se" s 



VI + 2x 2 s 2 



ds . 



(Bl) 
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Thus, Eq. (41a) holds for N = 0. The last expression in 
Eq. (Bl) is the same integral which appears in Eq. (26) of 
Hohenberg and Werthamer. 9 We next rewrite Eq. (37a) 
with respect to r\ jv in Eq. (40a) as 



(-% + I)/* 2 



(N = l) 



Vn 



(-TIN-! + VN^lm-2)/VNx 2 (N > 2) 



(B2) 

Using Eqs. (Bl) and (B2), it is easy to see that Eq. (41a) 
holds for N = 1. Proceeding to the general case, we as- 
sume that Eq. (41a) is valid for N < M — 1. We also 
remember the following properties of the Hermite poly- 
nomials: 

H N {s) - 2sH N _ 1 {s) + 2(N-i)H N - 2 (s) = , (B3a) 



Thus, we have established the first expression in Eq. 
(41a). The proof for the second expression proceeds in 
the same way by using partial integrations for rjo and 
?7at-2 in Eq. (B2). Equation Eq. (41b) can be proved 
similarly by induction, starting from Ri = 1 and using 
Eq. (37b). ' 



p CO 

/ s k H N (s)e- s2 ds = (k<N-l). (B3b) 
Jo 

is obtained explicitly by using Eq. (B2) as 



Then ijm is 
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